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^^esearch  conducted  under  this  contract  included  two  distinct 
projects.  The  first  was  concerned  with  using  simple  truss  models  to 
illustrate  certain  critical  aspects  of  nonlinear  material  behavior.  Results 
obtained  were  issued  in  two  reports  listed  as  items  1  and  ^  in  the 
bibliography  on  page  2.  This  project  was  carried  out  by  the  iVincipal 
Investigator.  Philip  6.  Hodge.  Jr. 


r 


The  other  project  was  concerned  with  using  a  constant-strain  finite 
element  method  to  find  complete  solutions  for  the  stress  and 
displacement  in  a  plate  with  a  centered  crack^his  research  was  carried 
out  by  James  Malone.  Research  Assistant.*  under  the  joint  direction  of  the 
Principal  Investigator  and  of  Professor  Robert  Plunkett,  ^he  elastic 
solution  was  obtained  firs^t  was  presented  in  the  report  listed  as  item 


3  in  the  bibliography.^e  elastic/plastic  solution  was  completed  after  the 
expiration  of  this  contract.  A  description  of  this  phase  constitutes  the 
main  body  of  this  Final  Report;  it  is  planned  to  submit  the  material  for 
publication. 
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rangements  are  also  considered.  Features  or  the  finite  element  algorithm 
e  tra  total  computing  time  are  discussed.  The  finite  element  program  is 
h  the  Cray- 1  computer  and  the  effect  of  vectorization  on  computational 


ABSTRAa 

In  this  paper,  the  finite  element  method  is  ^lied  to  a  center-cracked  plate 
subject  to  op^ira  mode  tensile  loading.  A  compfete  elastic-plastic  plane-stress 
solution  for  strain  hardening  materials  obeying  a  Von  Mises  yield  condition  and 
PrandtI-Reuss  stress-strain  relations  is  obtained  using  only  constant  strain  elements. 
An  accurate  representation  of  the  stress-strain  field,  even  at  distances  very  close  to 
the  aack  tip.  is  achieved  by  the  use  of  a  mesh  arran^ent  in  which  the  size  of  the 
elements  deaeases  in  a  geometric  series  as  the  aack  tip  is  approached.  The 
numerical  solution  is  compared  with  and  used  to  discuss  the  range  of  validity  of  the 
well  known  HRR  (Hutchinson-Rice-Rosengren)  aack  tip  solution  valid  fa  small  scale 

M.  The  influence  of  diffaent  amounts  of  hardenira  and  the  effect  of  changes  in 
h  arrangements  are  aiso  considered.  Features  or  the  finite  element  algorithm 
which  reduce  the  total  computing  time  are  discussed.  The  finite  element  program  is 
executed  on  the  Cray-t  computer  and  the  effect  of  vectaization  on  computational 
speed  is  discussed  fa  this  problem. 

1.  INTRODUCTION. 

In  a  recent  paper  [1]  we  have  shown,  fa  the  paticula  problem  of  a 
center-aacked  elastic  plate  unda  tensile  loading,  that  a  finite  element  famulation 
which  uses  only  constant-strain  elements  can  provide  an  accurate  representation  of 
the  stress-strain  field  even  in  the  neighbahood  of  the  stress  and  strain  singulaities 
at  the  aack-tip.  This  was  achieved  py  using  a  lage  number  of  triangula  elements 
arranged  in  a  mesh  in  which  the  size  or  the  elements  deaeases  in  a  geimetric  series 
as  the  aack  tip  is  approached.  The  present  paper  extends  this  tecnnique  to  handle 
nonlinea  material  behavior  and  discusses  the  resulting  complete  elastic-plastic 
solution  fa  aack  problems. 

Methods  commonly  used  in  nonlinea  finite  element  analysis  ae  discussed  by 
Bathe  [2]  and  Zienkiewicz  [3];  practical  procedures  fa  vaioys  plications  have  been 
presented  by  Bathe  and  Cimento  (4I  and  also  by  Bagan  et  al  [51.  In  an  elastic-plastic 
finite  element  analysis,  the  most  effective  way  of  dealing  with  the  mataial 
nonlineaity  is  by  an  inaemental  approach  in  whichihe  load  is  applied  in  a  number  of 
small  inaements.  within  each  inaement.  the  true  stress-strain  relations  are 
approximated  and  an  iterative  process  is  applied  to  ensure  that  the  equilibrium 
equations  are  satisfied  to  within  some  specified  level  of  accuracy.  At  each  itaative 
step  the  solution  of  a  banded  system  of  linea  algebraic  equations  is  required  which 
accounts  fa  a  lage  pation  of  the  computational  cost.  Fa  a  given  finite  element 
mesh,  the  accuracy  or  the  solution  can  be  improved  by  reducing  the  size  of  the  load 
inaements  and/a  imposing  a  tighta  convagence  aiterion.  Howeva.  this  may  cause 
a  substantial  inaease  in  the  total  computing  time  so  that  in  practice  a  balance  has  to 
be  struck  between  accuracy  and  computational  cost. 


Access  to  a  CRAY-1  supercomputer  has  made  it  feasible  for  us  to  consider  a 
technique  for  nonlinear  finite  element  anaiusis  which  uses  a  large  number  of  unknowns 
and  many  small  load  increments.  This  is  possible  because  the  computational  spe^  and 
memory  capacity  of  this  machine  are  much  greater  than  previous  generations  of 
computers.  For  example,  rates  of  69  MFLOPS  fmillion  floating  point  operations  per 
second)  have  been  recorded  [6.7]  on  the  C^Y-lS  for  the  solution  of  dense  systems  of 
linear  equations  of  order  300  usind  LU  decomposition  with  pivoting  ( without  resorting 
to  the  use  of  assembly  language),  mis  is  about  28  times  faster  than  the  IBM  3033  ana 
600  times  faster  then  the  w  11/760.  Rates  in  excess  of  HO  MFLOPS  have  been 
achieved  by  the  use  of  assembly  language  (8). 

In  this  paper,  we  consider  the  particular  problem  of  a  center-cracked  plate  subject 
to  mode-1  (tensile  opening  mode)  loading.  Small  strains  and  displacements  are 
assumed.  This  means  that  material  nonlinearities  but  not  geometric  nonlinearities  are 
considered.  This  problem  has  been  the  subject  of  considerable  interest,  in  the 
following  we  will  mention  only  some  of  the  many  analytic  and  numerical  papers  which 
have  applied  in  the  literature.  A  more  detailed  account  is  contained  in  the  review 
article  by  Rice  [9]. 

The  asymptotic  form  of  the  plastic  aack  tip  stress-strain  field  in  plane  strain  and 
plane  stress  has  been  established  analytically  by  Rice  [101.  Rice  and  Rosengren  [It], 
and  Hutchinson  [12.  13]  for  strain  hardening  and  perfectly  plastic  materials.  This 
soiutioa  frequently  referred  to  as  the  HRR  solution,  is  based  on  the  deformation  theory 
of  plasticity  and  is  valid  under  conditions  of  small  scale  yielding  [Hj.  It  will  be 
discussed  in  more  detail  later  in  this  paper. 

For  power  law  hardening  materials.  Tracey  [151  used  the  finite  element  method  to 
determine  the  plane-strain  ^ess  state  at  the  tip  of  a  crack  under  conditions  of  small 
scale  yielding.  He  used  special  singularity  elements  at  the  crack  tip  in  which  the 
displacement  shape  functions  were  cnosen  to  represent  the  form  of  the  HRR  solution; 
elsewhere  q-noded  isoparametric  elements  were  used.  Hilton  and  Hutchinson  [161. 
carried  out  a  plane-stress  finite  element  analysis  for  both  small  and  large  scale 
plastic  yielding.  Their  method  used  constant  strain  triangular  elements  togetn^  with 
a  special  singular  element  surrounding  the  aack  tip  in  which  the  HRR  solution  was 
embedded. 

We  have  recently  shown  [1].  fa  an  elastic  analysis,  that  constant  strain  elements 
can  be  used  to  obtain  an  acocurate  representation  of  the  stress-strain  field  in  the 
region  near  the  aack  tip.  Fa  the  elastic-plastic  analysis,  our.  ai^oach  has  also  been 
to  use  a  large  numba  of  constant  strain  triangula  elements  and  a  mesh  which  is 
simila  to  that  used  in  the  elastic  finite  element  analysis  [11.  Our  nonlinea  finite 
element  procedure  yields  a  complete  elastic-plastic  solution  at  a  reasonable 
computational  cost  when  executed  on  the  CRAY-1.  We  demonstrate  that  on  a  scale  in 
which  the  aack  length  equals  one  it  is  possible  to  ob^in  a  sufficiently  accurate 
solution  fa  the  stress  field  at  distances  as  close  as  10  from  the  aack  tip.  This 
solution  is  used  to  discuss  both  the  spatial  range  ova  which  the  HRR  solution  holds 
and  the  level  of  applied  load  at  which  the  HRR  solution  no  lor^  accurately  re^esents 
the  crack  tip  stress-strain  field.  Our  solution  provides  a  good  representation  of  the 
stress-strain  field  at  all  levels  of  applied  load. 

In  section  2.  we  desaibe  the  element  mesh  and  the  algaithm  which  implements 
our  elastic-plastic  finite  element  analysis.  The  numaicar  results  ae  discussed  in 
Section  3.  The  effect  of  diffaent  element  meshes  on  the  computed  solution  and  the 
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I  influence  of  hardening  are  also  reported  in  that  section.  In  Section  4,  computing  times 

for  different  meshes  and  for  different  amounts  of  hardening  are  discussed.  Features 
of  the  algorithm  including  vectorization  which  reduce  the  computational  cost  are  also 
desoribeo.  Finally,  some  conclusions  are  presented  in  Section  5. 

i  2.  METHOD  OF  SOLUTION. 

In  this  paper,  the  particular  problem  of  a  re^angular  plate  containing  a  centrally 
^  located  crack  subject  to  mode-f  (opening  mode)  loading  is  considered.  The  aack 

;  length  is  2a  and  the  plate  is  of  height  2n  and  width  ^a  (see  Fig.  I).  The  boundary 

I  conditions  along  the  top  and  bottom  edges  are  uniform  dispiacements  in  the  vertical 

I  direction  and  zero  traction  components  in  the  horizontal  direction.  Zero  tractions  are 

prescribed  on  the  remaining  edges  of  the  plate  and  along  the  crack  faces  SG.  By  the 
usual  symmetry  argument,  the  problem  can  be  reduced  to  that  of  solving  for  one 
quadrani  ABEF  (shac^  area)  of  the  plate  with  boundary  conditions  as  shown  In  Fig.  I . 

(a)  Finite  element  mesh. 

f 

]  The  finite  element  program  is  based  on  the  well  known  displacement  method  and 

uses  only  constant  strain  elements.  The  rectangular  domain  ABEF  (Fig.  I)  is 
disaetized  into  triangular  elements  by  a  mesh  generating  program  which  nas  been 
developed  so  that  different  mesh  arr^vements  can  be  automatically  produced.  A 
typical  mesh  arrangement  using  216  elements  is  shown  in  Fig.  2.  Meshes  ranging  in 
I  size  from  564  to  2064  elements  and  having  elements  as  small  as  10  a  at  the 

N  crack  tip  have  been  used  in  this  study. 


The  mesh  over  portion  ABCO  of  the  plate  is  formed  from  the  quadrilaterals  defined 
by  a  set  of  rectangular  rings  intersected  by  rays.  Let  M  numbers  r,  be  defined  by 


r„  =  a.  r,  =  o<r,*, .  t=M-I.M-2. ...  2.1  (2.  I) 

where  «  is  a  constant.  The  rings  are  a  set  of  nested  r«  x  Zr.  rectangles  with  lower 
corners  along  AGB  (see  Fig.  2)  at  distances  ±  r,  from  the  aack  tip  6.  Let  N-1 
equidistant  points  be  inserted  along  BC  and  the  same  spacing  continued  along  CD  and 
DA.  The  rays  are  strai^t  lines  from  G  through  these  points.  The  innermost  ring  is 
divided  into  4N^  equal  fscosceles  triangles  and  the  quadrilaterals  in  the  remaining 
rings  are  each  split  into  two  triangles  bu  their  diagonals  as  shown  in  Fig.  2.  Finally 
the  mesh  is  completed  by  filling  the  remaining  pation  DCEF  with  approximately  square 
rectmles  of  constant  size  which  are  split  into  triangles.  The  mesh  shown  in  Fig.  2  is 
fa  M=5  and  N=3. 

The  aspect  ratio  F  of  a  genaic  triangle  JKL  in  ring  \*  1  along  GB  is  given  by 
F=  |JK|/|JL|  =(r,^,-r,)/(r,/N)  =  Nll/oe-ll  (2.2) 

Hence  the  geometric  coefficient  oc  in  Eq  (2.2)  is  related  to  a  typical  aspect  ratio  by 


(2.S) 


We  have  found  that  good  results  ae  obtained  by  taking  F=l. 

A  similar  mesh  has  been  used  (II  in  an  elastic  study  of  an  infinite  plate  containing 
a  centrally  located  aack.  The  main  features  of  this  mesh  arrangement  are  mae  fully 
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discussed  in  Ref.  [I ). 

(b)  Material  properties. 

The  plate  material  is  assumed  to  be  homogeneous  and  initially  isotropic,  obeying 
Hooke's  law  in  the  elastic  range  with  initial  yielding  determineo  by  the  Von  Mfses 
yield  conditioa  After  yielding,  plastic  strain  increments  are  defined  by  an  associated 
flow  rule  with  linear  hardenira.  Isotropic  growth  of  the  yield  condit^  is  assumed. 

n 

apply.' '  These  assumptions'*  lead  to  the  well  known  PrandtI-Reuss  stress-stfain 
relaiions  in  the  incremental  theory  of  plasticity.  The  solution  presented  in  this  paper 
is  for  a  state  of  plane-stress.  It  is  also  assumed  that  both  strains  and  rotations  are 
small.  The  validity  of  these  assumptions  will  be  discussed  in  Section  3. 

The  elastic-plastic  stress-strain  matrix  0(<7)  which  relates  changes  in  stress  to 
changes  in  strain  at  points  of  the  material  which  have  uielded  and  are  loading  has  been 
derived  under  the  above  assumptions  by  Yamada  et  al  n?].  During  any  time  increment 
At  the  stress  increments  Ac  and  strain  increments  A£  are  related  by 


AC  =  D((y)  Ac 


(2.  A) 


where  AO  =  lAOv  .  AOu  ,  Ar  vn  F  .  Af  =  lAew  ,  Aft.  ,  Atf  wn  .  and  is  the 
engineering  shear  component  of  ^raia  If  the  materiaPis  eiasifc.  Over)  is  the  constant 
mArix 


I  0  0 

D(a)  =  [)c«E/(l-o2)  0  I  0 

^  0  0  (I-o)/2 


(2.  5) 


whether  or  not  earlier  plastic  behavior  has  taken  place.  If  It  is  plastic. 


D(a)  =  D^.(E/Q)  f 


(2.6) 


(On  ♦’b0wV(l*o)  R/2(l%)^2H/gCXl-o)  02 


-xy\«y 


where  a  is  the  deviatoric  stress  and  7  is  the  equivalent  stress. 


g  =  [(3/2)(0'||0’|j)l'^ 

P  =  (2H/9E)g2*r^2/(l*o) 
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4.  IMO 


Q  »  R  ♦  2(1-d*)P 
R  *  0„**  ♦2T)Og*Cy‘»  0y’2 
H  =  EEt/(E-Et) 


(2.  7) 


where  o  is  the  Poisson's  ratio.  E  is  Young's  Modulus,  and  Ij  is  the  plastic  modulus. 

A  plot  of  equivalent  stress  7  versus  equivalent  strain  r  for  a  material 
Sr^d  assumptions  is  shown  in  Fig.  3.  The  equivalent  strain  e  is 

r  =  I(3/2Xf’jjf‘|j)l‘^  (2.8) 

where  £'jj  are  the  deviatoric  components  of  strain.  In  practice  a  piecewise  linear 
stress-strain  law  can  be  determined  from  data  obtained  by  a  simple  tension  test  for  a 


where  £jj  are  the  deviatoric  components  of  strain.  In  practice  a  piecewise  linear 
stress-strain  law  can  be  determined  from  data  obtained  by  a  simple  tension  test  for  a 
real  material.  We  have  chosen  the  most  simple  case,  i.e.  a  bilinear  stress-strain  law. 
but  the  above  discussion  can  be  modified  in  an  obvious  manner  to  handle  any  piecewise 
linear  curve. 


(c)  Nonlinear  finite  element  analysis. 

Since  the  strain  is  constant  in  each  element,  the  column  matrix  of  element  strains 
£  and  the  column  matrix  of  nodal  displacements  u  are  related  by 


£  s  B  u 


(2.9) 


where  the  constant  matrix  6  depends  only  on  geometry  and  not  on  material  behavior. 

The  nodal  force  matrix  F  and  the  stress  matrix  a  are  related  by  the  similar 
relation 

F  =  Ao  (2.10) 

where  A  also  depends  only  on  geometry.  Equations  (2.  9)  and  (2.  10)  must  be  valid  for 
either  total  or  increment^  benavior. 


During  any  incremental  time  step  At,  the  stress  and  strain  increments  are  related 


a<j  =  0(0)  * ,  [)((,.) 


(2.  II) 


Where  0**  is  some  mean  value  of  the  stress  state  during  the  strain  interval  A£ . 


K(0)Au  =  AF 

where  the  'stiffness*  matrix  K  is  given  by 

K  (0)  =  A  D(0)  B 


(2.  12) 


(2.  13) 


Due  to  the  incremental  nature  of  the  stress-strain  relations  (see  Eq.  (2.  1 1)).  this 
problem  is  best  solved  by  an  incremental  procedure  in  which  the  load  is  applied  in  a 


series  of  increments  .  Thus  the  load  at  the  end  of  increment  M  is 

F„=F„.,  M=!.2 . q  (2.  M) 

where  Fq  .  F^  are  the  initiai  and  final  loads  respectively. 

The  problem  can  be  posed follows:  Given  a  complete  solution  u .  €  m.i  . 
a  if.,  at  the  end  of  the  (n-l)^  step  corre^oiYling  to  the  apt^ed  extemarioad  F>!.| . 
fiiVl  the  complete  solution  u  M  .<h  at  the  end  of  the  fP  step  corresponding  to 
the  external  load  F,^ .  By  a  complete  Solution  .  we  mean  column  matrices  for  nodal 
displacements,  .nodal  forces,  element  strains,  and  element  stresses.  The  process 
during  the  n  step  is  to  obtain  incremental  changes  Au  .  Ac  ..  .  Aa  »  in 
displacements,  strains,  and  stresses,  respectively,  corresponding  to  the  increment  in 
load  which  satisfy  Eqs.  (2.  9)  to  (2.  1 1).  i.e. 


Ac  ff  =  BAu  ff 
Aa„sD(o„  "  )  Ac„ 


(2.  15) 
(2.  16) 
(2.  17) 


where  a„  ~  0^.1  AOm  gIm  **  is  some  mean  value  of  the  stress  state  during  the 
strain  inwvai  Ac  w^l  be  defined  later). 

The  first  step  is  fully  elastic  and  the  solution  is  obtained  by  direct  solution  of  the 
elastic  finite  element  equations.  For  each  later  step  an  iterative  method  is  used  to 
solve  the  kinematic,  constitutive,  and  equilibrium  equations  (2.  15,  -16.  -17)  at  step 
M.  The  first  estimate  Au  n'*'  is  given  by 

Au„^^^=  AUrt_|  (2.  18) 

where  At  n-i  is  the  change  in  displacement  which  occurred  during  the  previous  (M-l) 
step.  Next  Ac  and  AcJ”  are  obtai^d  in  order  from  Eqs.  (2.  15)  and  (2.  16) 
using  a„.,  as  the  lnitial,roean  stress  ;  the  first  total  stress  estimate  aJ” 

is  o^ainM  by  adding  to  the  initial  value  a^.,. 

In  general  oJ'l  will  not  satisfy  equilibrium.  We  begin  the  iterative  process  by 
defining  a  residua)  force: 

P„''>=F„ (2.19) 
and  obtain  a  correction  to  the  displacement  increment  field  by  solving 

ic„$u  „<•>=?„">  (2.20) 


Where 


Kfi  ®  A  D(o^f )  B 


(2.21) 


adding  Su  to  Au  and  Ac  is  obtlined  from  (2. 15).  A  mean  stress  a 


is  defined  by 


(2.  22) 


a„*(2)=(a„.,*a„o))/2 

and  is  then  given  by  (2.  16). 

This  iterative  process  is  repeated  untii  the  residuai  forces  at  the  i^^ 
iteration  satisfy 


1 1  P„«>  1 1  <  Ta  (2.  23) 

wh^e  1 1  PJ')  1 1  is  the  magnitude  of  the  largest  component  (in  absolute  value)  of 
Pm”^  over  all  the  nodes  and  TOL  is  som^eset  convergence  tolerance. 

The  iterative  procedure  at  the  step  can  be  summarized  in  the  following 
algorithm: 


Au  =  Au  ♦  Su 
A£  Au„«*» 

♦<?„«>  )/2 

=  ♦DC  <?„“«♦>>)  A€ 

P„<‘-«)=F„-Acr„«-n 


ic„su„«*i)=P„(‘-’) 


where  8u  - =  0  ,  Au  Au ,0  -W)  =  ^  ^  o,l,2,3 . It  can  be 

seen  that  inis  approach  is  a  modif  leo  form  of  the  Newton-Raphson  method  for  solving  a 
non>l inear  system  of  equations  (see  Refs.  12,  3]). 


(d)  Load  step  size. 


One  of  the  advantages  of  using  constant  strain  elements  is  that  at  any  given 
applied  load  the  stress  tnroughout  each  element  will  be  constant.  Further  Tor  the 
problem  considered  here,  no  unloading  of  yielded  elements  occurs.  Therefore,  at  any 
given  load  each  element  is  either  elastic  and  has  never  yielded  or  is  plastic  ana 
loading.  Our  numerical  procedure  is  to  terminate  a  load  step  when  ^  elastic  element 
reaches  yield.  When  this  procedure  is  used  each  element  will  remain  either  elastic  or 
plastic  tnroughout  the  entire  load  step. 


Since  any  single  iteration  is  a  strictly  linear  process,  the  load  Ff.  at 
s^gp  js  terminated  can  be  easily  estimated  in  the  following  way. 
beginning  of  the  M  load  step  we  apply  some  arbitrary  load  increment 


_  3f  the  n  load  step  we  apply  s^e  arbitrary  loa< 
proceed  to  compute  Au  J” ,  A«  ,  Aoh'”  in  the  manner 
for  each  elastic  element  p  we  compute  a  ^lar  factor 


which  the 
At  the 
AFf.  and 
described  above.  Then 
which  satisfies 


(2 . 25) 


where  Y,  is  the  initial  yield  stress.  Based  on  the  stress  change  AoJ'*  we 
predict  that  the  element  «  tor  which 
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(2.26) 

will  be  the  element  most  likely  to  first  reach  yield  duhng  the  current  load  step.  The 
entire  incremental  solution  can  now  be  scaled  by  the  factor  7f„ 

(2.27) 

so  that  element  «  will  have  Just  reached  yield  when  the  load  is  incremented  by  an 
amount 


AF„  AF„  (2.28) 

At  each  iteration,  improved  estimates  of  the  load  increment  AF-  requir^  to  have 
element  a  just  reach  yield  can  be  obtained  by  computing  scaling  fxtprs  for 
element  «  corresponding  to  the  stress  change  Aa^  «  - 

It  was  found  that  computing  scaling  factors  for  only  the  first  three  iterations  in 
each  load  step  was  the  most  efricimt  procedure.  The  l(Sd  increment  AFg  is  then  held 
fixed  for  the  remaining  iterations  of  the  load  step.  At  the  end  of  the  load  step,  t^^  yield 
condition  may  not  be  satisfied  exactly  by  element  «  but  this  is  taken  care  of  by  the 
use  of  a  “smeared”  yield  condition. 

The  idea  of  a  smeared  yield  condition  is  based  on  a  technique  used  by  Yamada  et  al 
(17}  and.  in  a  different  context.  1^  Hodge  and  Van  RIj  1 18.  19]  to  substantially  reduce 
the  total  number  of  load  steps.  At  the  end  of  each  load  step,  any  elastic  element  3 
for  which  the  equivalent  stress  9  satisfies 

0,99  Y,  <^<I.02Y,  (2.29) 

has  its  yield  stress  Y|  redefined  so  that  Y|  =  7 . 

This  means  that  elements  which  have  not  yet  reached  yield  but  which  are  close 
enough  to  satisfy  (2.  29)  can  be  treated  as  plastic  in  the  next  step.  This  avoids  the 
need  for  one  or  more  additional  steps  to  brira  these  elements  to  yield.  It  also  means 
that  the  solution  at  the  end  of  the  step  can  be  accepted  even  if  some  elements  exceed 
the  yield  stress,  provided  (2.  29)  is  satisfied.  The  alternative  would  be  to  repeat  the 
load  step  using  an  improved  estimate  for  the  initial  size  of  the  load  increment. 

The  main  features  of  the  program  have  been  described  above.  A  more  detailed 
desaiption  of  the  program  algoriCnm  may  be  found  in  Ref  120}  which  contains  the 
program  documentation  and  Fortran  code. 


3.  RESULTS  AND  DISCUSSION. 

In  this  section  an  assessment  of  the  accuracy  of  the  solution  particularly  in  the 
region  close  to  the  crack  tip  will  be  given.  The  effect  of  different  mesh  arrangements 
on  the  solution  will  also  be  discussed.  Comparisons  will  be  m^  between 
elastic-plastic  solutions  corresponding  to  different  amounts  of  hardening  and  also 
with  the  purely  elastic  solution.  The  numerical  elastic-plastic  solution  in  the  region 
surrounding  the  crack  tip  will  be  compared  with  an  analytic  solution  which  is  valid 
under  conditions  of  small  scale  yielding  in  which  the  plastic  zone  is  small  compared  to 
the  plate  dimensions. 


(a)  Finite  element  solutioa 

.  '*'0^09'*  Modulus  E  =  185.1  G  Pa,  Poisson's  ratio  o  =  0.3  . 

IniUal  yield  stre»  Y  =  225.-4  M  and  ob^s  a  bilinear  stress-strain  law  with  plastic 
mo(ajlus  E  j  =2.372  6  Pa.  The  plate  dimensions  are  h=250  mm  and  a=^  mm  (see  Fig. 
I).  This  pvticular  oeometru  and  material  properties  were  selected  so  that  the 
n^vical  solution  could  also  m  compared  with  the  experimental  results  of  Yagawa  et 
ai.  u  I  ]. 

They  reported  the  increase  of  a  gage  length  in  a  cmtrally  aacked  thin  (5  mm  thick) 
aluminum  plate  as  a  function  of  load.  The  gage  points  straddling  the  crack  were 
located  at  distances  of  80  mm  directly  above  and  below  the  center  of  the  plate.  Our 
finite  element  solution  was  obtained  using  mesh  A  (see  Table  I)  and  agrees  closelu 
with  the  experiment  as  shown  in  Figure  4.  A  dimensionless  load  f  has  been  defined  by 
dividing  by  the  Initial  yield  load  of  the  uncracked  plates  f*  p/A.Y  where  A,  is  the  area 
of  ^  lop  edge  of  the  plate.  Henceforth.  *  applied  load  *  reisers  to  f.  Ya^a  et  al 
(21)  found  a  gage  point  displacement  of  4.43  mm  at  an  applied  load  of  0.87  at  which 
loM  the  crack  began  to  grow,  whereas  for  the  same  displacement  our  numerical 
solution  predicts  an  applied  load  of  0.93,  which  is  7%  greater  than  the  measured  value. 

The  plastic  zones  at  different  levels  of  applied  load  are  shown  in  Figs.  5  and  6.  It  can 
be  seen  that  the  overall  shape  of  the  plastic  zone  is  somewhat  influenced  by  the  level 
of  applied  load.  Our  numerical  results  show  that  the  plastic  zone  has  a  radius  of 
about  p. la  in  the  region  ahead  of  the  aack  tip  at  an  ^lied  load  of  0.3  .  The  plastic 


and  extends  to  cover  60x  of  this  area  at  a  load  of  0.9 . 


Tlw  equivalent  stress  9  computed  at  an  applied  load  of  0.3  along  5  rays  radiating 
from  the  aack  tip  is  plotted  versus  Ira  (r/a)  In  Fig.  7  where  r  and  8  are  pola 
coadinates  referred  to  the  aack  tip.  The  radial  distance  r  is  measured  from  the 
tip  to  the  centroid  of  the  element  in  which  ^  is  computed.  .  Each  curve 
displays  a  fairlu  abrupt  kink  indicated  fa  example  by  the  arrow  in  Fig.  7  fa  the  ray 
0  =168  ®  .  Along  each  ray  the  kink  ococurs  at  the  elastic-plastic  boundary.  The 
stresses  drop  off  sharply  ova  a  shat  distance  beyond  the  elastic-plastic  bounday. 
Simila  behavia  is  observed  at  otha  levels  of  applied  load. 

.  applied  load  Inaeases  above  0.5  and  the  plastic  zone  extends  aaoss  the 

plate,  it  is  found  that  the  equivalent  stresses  in  the  region  close  to  the  aack  tip 
(r/a<IO  "2 )  are  greata  than  but  approximately  propationai  to  the  equivalent  stresses 
obtained  from  the  purely  elastic  solutioa  The  f^ors  of  propationality  ae  about  2.6 
and  1.4  at  applied  loara  of  0.6  and  0.9,  respectively.  This  means  that,  at  applied 
loads  greata  than  0.5  the  stress  field  dispiras  a  1//r  singulaity  at  the  aack  tip. 
The  detailed  behavia  of  the  stress  field  at  applied  loads  below  f=  0.5  will  be 
discussed  lata  in  subsection  (d). 


®  shows  the  ^  components  of  strain  on  a  log-log  scale.  The  strains  along 
the  49  ®  r^  show,  a  l//r  variation  at  the  aack  tip  ova  the  ranra  r/a  <  2.0  x  10 
.  Along  this  ray.  the  elastic-plastic  boundary  is  located  at  r/a  =  0.1 5.  Ova  the  range 
r/a  >  0.15  the  ^  strain  components  ae  greata  than  those  which  would  be  obtained 
from  a  purely  erastlc  solutioa  Fa  example,  the  percentage  diffaence  between  the 
strains  from  the  elastic-plastic  and  elastic  solutions  dkreases  from  10X  just 


outside  Uw  elastic-plastic  boundary  to  less  than  3X  near  the  edge  of  the  plate  (for  r/a 
>  0.9) .  As  the  applied  load  is  increased  the  distance  from  the  aadc  tip  over  which  e. 
dispiMs  a  1/  /r  variation  increases  in  extent  up  to  a  maximum  distance  of  about 
r/a  =  0.02  along  the  49  ^  ray  at  an  allied  load  or  0.9.  The  behaviour  of  the  strain 
along  other  rays  is  qualitatively  similar:  the  largest  c,  strain  components  occur 
alorig  the  ray  for  which  0=61*  ^ 

We  have  also  considered  a  material  with  a  higher  amount  of  hardening  (  E  r  = 
18.51  6  Pa  ).  The  behavior  of  the  stress-strain  field  for  this  material  is 
qualitatively  similar  to  that  obtained  for  the  lower  hardening  material  (E  ^  *  2.372  G 
Pa)  which  has  been  discussed  above. 

(b)  Effect  of  different  mesh  arrangements. 

The  presence  of  a  singularity  in  the  stress-strain  field  at  the  aack  tip  requires 
that  the  arran^ent  of  elements  in  the  mesh  must  be  carefully  selected  if  an 
accurate  solution  is  to  be  achieved  at  a  reasonable  computational  cost.  The  elastic 
finite  element  solution  for  a  crack  problem,  using  a  mesh  similar  to  that  described  in 
Sec.  2  is  discussed  in  Ref.  [1].  It  was  found  that  increasing  the  number  of  rings  in  the 
mesh  produced  more  accurate  results  in  the  region  close  to  the  crack  tip.  whereas  an 
increase  in  the  number  of  rays  gave  more  accuracy  over  the  rest  of  the  plate.  It  was 
also  shown  that  for  a  given  mesh,  increasing  the  number  of  rings  while  holding  the 
number  of  rays  fixed  ( thereby  increasing  the  density  of  elements  only  at  the  crad  tip 
)  produced  less  than  O.IX  improvement  in  the  accur^  of  the  solution  over  the  range  r 
>l00r.  where  ri  is  the  position  of  the  innermost  ring  in  the  original  mesh.  We  now 
show  that  similar  observations  hold  for  the  elastic-plastic  problem. 

Four  specific  mesh  arrangements  have  been  considered  (Table  I).  Solutions  for  the 
high-hardening  material  (Er  =16.51  G  Pa  )were  obtained  for  each  mesh  over  a  range  of 
loads  up  to  0.9  .  At  this  load  the  plastic  zone  extended  to  touch  almost  all  of  the  edge 
EF  (see  Fig.  2).  To  facilitate  a  comparison  with  the  results  of  Ref.  [II  .  we  will 
compare  the  vertical  components  of  displacement  v  at  nodes  along  the  aack  face 
obtained  from  each  mesh.  Similar  behavia  is  observed  fa  v  aiong  other  radial 
directions  and  also  fa  the  strain  field. 

Meshes  A  and  B  v  ^  iiosen  to  study_the  effect  of  changing  the  number  of  rings 
while  keeping  the  numbers  of  rays  fixed.  The  solutions  obtained  from  meshes  A  and  B 
show  no  oiffaence  in  v  over  the  range  10  <  r/a  <  1  at  all  levels  of  loading. 

However,  over  the  range  10  <  r/a  <  10  the  vertical  displacements  v  obtained 

from  mesh  B  are  less  than  those  obtained  from  mesh  A.  The  diffaence  between  the 
two  solutions  inaeases  as  the  aack  tip  is  approached.  Fa  example,  at  an  applied 
load  of  0.3  the  differences  in  v  are  appproximately  0.5X,  5.5X,  and  IdX  at  r/a  equal 
to  10  .  10  .  and  3.33x10  ( the  position  of  tne  node  on  the  aack  face  adjacent 

to  the  aack  tip  in  mesh  B).  respectively. 

The  effect  of  inaeasing  the  density  of  the  elements  over  the  entire  plate  by 
doubling  the  number  of  rays  while  holdira  the  number  of  rings  fixed,  was  studied  by 
comparing  meshes  C  and  D.  The  vertical  displacements  v  along  the  aack  face  obtained 
from  the  coarse  mesh  C  are  smalla  than  those  obtained  from  the  refined  mesh  0.  Fa 
example,  at  an  applied  load  of  0.3  the  differences  in  v  were  approximately  I.4X,  2%, 
3X,  and  7X  at  r/a  equal  to  I.  10  •* .  10  ^  .  and  4.6x10  .  respectively. 

An  assessment  of  the  accuracy  of  the  finite  element  solution  fa  the  elastic 
problem  [11  was  made  possible  by  a  comparison  with  an  exact  analytic  solution.  It 


-  v 


was  found  that  the  elastic  strain  and  displacement  fields  obtained  by  the  use  of  mesh 
A  were  in  very  good  agreement  with  the  analytic  solutioa  For  example,  the  error  in 
vertical  di^lacement  v  along  the  aack  face  was  2X,  4X,  6X,  and  i4X  at  r/a  equal  to 
1.  10^.  10^.  and  lO"*.  resp^lvely.  ^ 

In  the  absence  of  a  complete  analytic  solution  for  the  elastic-plastic  problem,  it 
is  not  possible  to  determine  the  accuracy  of  the  finite  element  solution  directly,  as 
was  done  for  the  elastic  problem.  However,  the  observations  which  have  been  m^  on 
the  effect  of  different  mesh  arrargements  for  the  elastic-plastic  problem  are 

Salitativelu  the  same  and  in  close  quantitative  agreement  with  those  made  in  Refill 
r  the  elastic  problem.  Given  this  agreement,  we  can  infer  from  the  results  of  [l] 
that  mesh  A  should  also  be  ejected  to  provide  an  element  arrangement  for  which  an 
accurate  elastic-plastic  solution  can  be  obtained  over  the  range  r/a  >  10  .  This 

assertion  Is  borne  out  later  in  the  paper  when  the  numerical  solution  is  compared  with 
the  analytical  HRR  crack  tip  solution  under  conditions  of  small  scale  yielding. 

For  the  elastic-plastic  problem,  it  is  found  (see  Section  5)  that  the  computational 
cost  increases  by  a  factor  of  about  10  when  the  number  of  rays  in  the  mesh  is  doubled 
(i.e.  ch^ing  N=3  to  N*6) .  However,  from  the  elastic  results  HI.  it  would  be  expected 
that  this  change  in  the  mesh  would  only  slightly  improve  the  accuracy  of  the  solutioa 
e.g.  by  approximately  li.  2X,  3X.  and  6%  in  v  at  r/a  equal  to  1.10  .  10  .  and  10 

.  respectively.  It  was  concluded  that  mesh  A  can  oe  expected  to  provide  the  best 
element  arrangement  in  terms  of  balancing  computational  cost  and  accuracy  over  the 
range  r/a  >10  ^  for  the  elastic-plastic  analysis. 

(c)  Small  scale  yielding  and  the  HRR  solutioa 

The  analytic  crack  tip  HRR  solution  (10.1 1.  and  121  is  based  on  a  deformation  theory 
of  plasticity  and  is  valid  under  conditions  of  small  scale  yielding.  The  term  small 
scale  yielding  refers  to  the  situation  in  which  the  applied  load  is  sufficiently  low  so 
that  the  size  of  the  plastic  zone  is  small  compared  to  the  length  of  the  crack:  it  is 
small  enough  that  the  plastic  zone  is  embedded  in  an  elastic  field  governed  by  the 
dominant  l7/r  term  in  the  asymptotic  elastic  series  solution.  In  obtaining  the  HRR 
solutioa  the  1//r  elastic  term  is  the  assumed  boundary  condition  for  large  r. 
However,  the  HRR  analysis  cannot  predict  how  large  the  plastic  zone  may  become  so 
that  the  t//r  elastic  term  is  still  a  good  appx’oximation  for  the  solution  in  the 
region  surrounding  the  plastic  zone.  Further,  the  HRR  solution  represents  the 
elastic-plastic  solution  only  over  a  small  region  of  the  elastic  zone  located  at  the 
crack  tip.  The  extent  of  this  region  cannot  be  determined  from  the  HRR  analysis. 

Our  numerical  results  provide  a  solution  over  the  entire  plastic  zone  and  are  based 
on  an  incremental  flow  theory  of  plasticity.  It  has  been  shown  [221  that  a  solution 
obtained  using  deformation  Theory  will  oe  similar  to  that  obtained  using  an 
incremental  flow  theory  of  plasticity  provided  the  condition  of  proportional  stressing 
is  satisfied.  It  has  been  pointed  out  (131  that  proportional  stressing  can  be  expected 
to  hold  under  the  assumption  of  small  scale  yielding  for  a  material  obeying  a  bilinear 
stress-strain  law. 


For  such  a  material  Hutchinson  [121  has  shown  that  the  radial  and  angular 
variation  of  the  stresses  in  the  HRR  solution  has  the  same  form  as  the  dominant  term 
in  the  elastic  solution,  in  the  HRR  solution  the  equivalent  stress  at  a  point  (r.  6) 
can  be  represented  by 

(cos2  (e/2)  ♦  (3/4)  Sin2  eV^  (3.  I ) 
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Where  r  and  6  are  polar  coordinates  referred  to  the  crack  tip.  The  amplitude  K  is 
given  by 


K  =  (Et/E)‘'2|c,  (3.2) 

where  K.  is  the  elastic  stress  intensity  factor  which  replaces  K  in  Eq.  (3.  1)  to  yield 
the  elastic  singular  term. 

Confidence  in  the  accuracy  of  the  numerical  solution  in  the  crack  tip  region  can 
be  gained  from  a  comparison  witti  the  HRR  solution.  On  the  other  hand,  the  accurate 
numerical  solution  can  be  used  to  assess  the  rar^e  over  which  the  Hfw  solution  Is 
valid  at  different  levels  of  applied  load  and  also  to  determine  the  behavior  of  the 
solution  over  the  remainder  or  the  plastic  zone  in  which  the  HRR  solution  does  not 
hold.  In  addition,  the  level  of  applied  load  up  to  which  the  HRR  solution  provides  an 
accurate  description  of  the  crack  lip  stress  field  can  be  estimated.  These  topics  will 
be  the  subject  of  the  discussion  that  follows. 

(d)  Comparison  with  the  HRR  solution  for  small  scale  yielding. 

The  elastic  stress  intensity  factor  K,  [used  to  compute  .  see  Eqs.  (3.  t)  and 
(3.  2)1  has  been  determined  from  the  elastic  finite  element^lution  by  the  method 
described  in  Refs.  [1  and  23).  It  has  been  shown  [1]  that  K,  computed  by  this  method 
will  be  accurate. to  within  3%  of  the  correct  value.  For  this  problem,  we  have 
computed  K,  //a  Y  *  1.23  .  We  compare  the  HRR  and  finite  element  solutions  for  the 

hardening  material  by  examinirg  the  ratio  X  *  ^  3long  the  49°  ray  where 

o  is  the  equivalent  stress  obtained  from  the  numerical  solufion.  The  results  are 
shown  in  Fig.  9. 

For  the  elastic  finite  element  analysis  we  know,  by  considering  the  analytic 
solution  for  ^  (given  by  the  1/  /r  elastic  term  In  the  region  near  the  crack  tm)  , 
that  the  error  in  o  increases  as  the  aack  tip  is  approached.  The  error  is  about  20% 
In  those  elements  nearest  the  crack  tip  but  is  less  than  15*  for  r/a  >  10"^  .  As 
discussed  in  subsection  (b)  above,  similar  behaviour  regarding  the  accuracy  of  9  can 
be  expected  from  the  elastic-plastic  finite  element  solution.  As  a  result,  we  can  be 
reasonably  confident  that  the  difference  between  7  and  at  points  for  which  X 
<  0.8  can  not  be  caused  solely  by  inaccuragj  in  the  numencal  solution.  The  HRR 
solution  provides  an  accurate  re^esentation  orthe  near  tip  stress  field  over  the  rarae 
r/a  <  p  where  p  is  unknowa  However,  from  Fig.  9  we  can  easily  determine  an  upper 
bound  on  p  by  using  the  criterion  tr»t  the  HRR  solution  can  be  considered  to  represent 
the  stress  field  oniy  at  points  for  which  X  >  0.6  . 

By  comparing  the  X  curves  for  f  equal  to  0.07, 0.22,  and  0.3  in  Fig.  9.  it  can  be 
seen  that  p  increases  in  magnitude  as  the  applied  load  is  increased.  However,  these 
curves  show  that  the  HRR  solution  represents  the  stress  field  only  over  a  very  small 
portim  of  the  plastic  zone.  For  example,  at  f=  0.3  the  elastic-plastic  boundary  along 
l^j49°  ray  Is  at  r/a  «  0.15  but  the  HRR  solution  is  certainly  not  valid  b^ono  r/^  = 

^  As  the  load  is  increased  above  0.45  ,  the  numerical  solution  for  the  stress  field 
at  the  aack  tip  begins  to  differ  inaeasinglu  from  the  HRR  solution.  Fa  example, 
consider  the  curve  caresponding  to  f  -  0.6.  This  is  expected  because  the  plastic  zone 
has  extended  to  the  edge  of  the  plate  (see  Fig.  5)  so  that  it  can  no  lon^  be  regarded 
as  embedded  In  an  elastic  field  ^wsrned  by  the  l//r  term.  Thaefae,  the  smalfscale 
yielding  conditions  assumed  fa  the  HRR  solution  are  no  lon^  true  at  these  higha 
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applied  loads.  From  our  numerical  results  we  conclude,  for  this  particular  aadc 
problem,  that  the  HRR  solution  does  not  represent  the  stress  field  anywhere  in  the 
vicinity  of  the  aack  tip  when  the  applied  load  f  exceeds  about  OS.  ft  can  be  seen 
from  tra  curve  for  f  =  0.6  that  the  stress  field  displays  the  l//r  singularity  over 
the  range  r/a  <  10'^  but  its  amplitude  is  no  longer  determined  by  K  in  Eq.(3. 2). 

Similar  behaviour  is  exhibited  by  the  elastic-plastic  solution  for  the  high  hardening 
material  (Er  =  18.51  6  Pa).  For  this  material.  It  is  found  that  the  range  over  which 
the  HRR  solution  represents  the  crack  tip  stress  field  is  slightly  larger  in  size  than 
that  found  for  the  low  hardening  material  at  the  same  applied  load. 

(e)  validity  of  the  small  strain  and  small  rotation  assumptions. 

The  components  of  strain  are  shown  in  Figure  8  for  the  low-  hardening  material 
at  an  applied  load  of  0.3  .  The  strains  along  the  ray  are  greater  than  0.1  over 
the  range  r/a  <10  .  it  is  clear  that  the  assumption  of  small  strains  is  violated  in 

the  region  near  the  aack  tip.  The  rotations  of  the  elements  in  this  region  are  also 
large. 

A  famulation  which  accounts  for  large  strains  and  large  rotations  would  be  a 
mae  appropriate  model.  However,  in  our  discussion  (subsection  (b))  on  the  effect  of 
the  different  meshes,  it  was  shown  that  changes  in  the  solution  over  the  region  near 
the  aack  tip  do  not  produce  significant  changes  in  the  region  away  from  the  aack  tip. 
The  large  strains  are  limited  to  a  very  small  region  at  the  aack  tip.  it  might  then  be 
expected  that  a  large  strain  formulation  would  not  significantly  alter  the  solution  over 
much  of  the  remainder  of  the  plate. 

The  same  conclusions  cannot  be  drawn  at  higher  levels  of  applied  load.  e.g.  f=  0.5 
or  greater,  for  which  the  plastic  zone  has  readied  the  outer  edge  of  the  plate.  At 
these  loads  the  region  of  large  strains  has  inaeased  in  size  (fa  example,  r/a  <  10 
at  f=  0.6)  to  the  extent  that  a  large  strain  solution  would  be  expected  to  differ  from 
the  small  strain  solution  ova  much  of  the  plate. 


4.  COMPUTING  TIME  AND  VEaORIZATION. 

In  our  finite  element  analysis,  the  use  of  a  large  numba  of  degrees  of  freedom, 
many  small  load  steps,  and  a  tight  convergence  aiiaion  is  feasibre  because  of  the 
speed  of  the  CRAY- 1  computa.  Fa  example,  a  typical  load  step  involves  the  solution 
of  1300  linear  equations  and  requires  about  30  itaations  fa  convagence:  the 
computations  fa  such  a  step  take  only  1.2  c  p.  u.  seconds. 

Fa  both  the  low  and  the  high  hardenira  mataials  the  size  of  the  load  inaements 
when  using  mesh  A  ranged  from  Af=  10  at  lower  applied  loads  up  to  Af=  10  at 
higher  loads.  The  tolaance  used  in  the  convagence  aitaion  (2.23)  was  set  at  TOL= 
l(f  and  was  held  fixed  fa  each  load  step.  At  this  tolaance  the  solution  fa  u  .  6 . 
a  and  nodal  faces  F  convaged  to  3  places  of  decimals. 

Information  about  the  computations  perfamed  using  diffaent  meshes  fa  both  the 
low  and  high  hadening  materials  is  displayed  in  Table  2.  Fa  the  same  mesh,  the  total 
c.  p.  u.  time  fa  the  low-hadening  mataial  is  about  twice  that  required  fa  the 
hign-hadening  mataial.  Fa  a  patTcula  mataial.  the  numba  of  itaations  required 
pa  step  is  independent  of  the  choice  of  mesh  and  is  governed  only  by  the  amount  of 
nonlineaity  present  in  the  problem.  i.e.  on  the  extent  of  the  plastic  zone  at  the 
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An  important  feature  of  the  CRAY- 1  is  its  vector  hardware  which  enables  the 
machine  to  execute  computational  processes,  in  which  the  same  arithmetic  operation 
is  performed  on  each  element  or  pair  of  elements  from  an  ordered  set.  by 
'vectorization*.  Vectorized  computations  are  performed  at  a  siignificant  increase  in 
sp^  compared  to  the  more  conventional  "scalar*  mode  in  which  computations  are 
pmormed  sequentially.  Guidelines  for  writing  FORTRAN  programs  which  make 
efficient  use  or  vectorization  can  be  found  in  Refs.  [24-26]. 

In  an  experiment,  vectorization  was  temporarily  turned  off  for  a  run  using  mesh  A 
so  that  vector  and  scalar  speeds  could  be  compared  for  our  code.  The  computations 
involved  in£qs.  (2.  15).  (2.  16).  (2.  17).  and  (2.  22)  were  carried  out  by  vectorization 
at  speeds  which  were  3.  5.  7.  and  10  times  heater,  respectively,  than  those  attained 
in  scalar  mode.  Assembly  of  the  global  stiffrass  matrix  and  load  column  matrix  is  an 
inherently  scalar  process,  the  small  amount  of  vectorization  which  could  be  achieved 
resulted  Tn  speeds  which  were  only  1.5  times  faster  than  scalar  mode. 

It  is  commcn,[2.31.  when  using  a  modified  Newton-Raphson  scheme,  to  obtain  an 
initial  guess  Au  for  the  change  in  displacement  during  the  h  ”  step  by  solving 
for  Au  from 

K„Au„t»=AF„  (4.1) 

However,  In  our  algorithm,  the  Initial  guess  Au  Is  given  by  the  displacement 
increment  Au  Trom  the  previous  step  (see  Eq.  2.  16).  we  found  that  this  simple 
dunm  substamially  reduced  the  number  of  iterations  required  for  convergence  in  each 
step.  As  a  result  the  total  c.  p.  u.  time  was  reduced  by  a  factor  of  about  2. 

In  the  modified  Newton-Raphson  Iterative  scheme  which  we  have  discussed  in 
Section  2,  the  stiffness  matrix  K  Is  updated  (using  the  current  state  of  stress)  at  the 
beginning  of  each  load  step.  An  alternative  Is  only  to  uDdate  K  periodically  and  to 
reCaIn  the  same  factored  form  (R^  R  Cholesky  algorithm)  for  a  number  of  load  steps. 
For  this  problem,  we  e^merimented  with  a  periodic  updating  strategy  based  on  the  rate 
of  convergent  of  the  iterative  process.  However,  we  found  that  this  strategy  was 
difficult  10  implement  successfully  in  our  program.  For  example,  our  best  attempts 
reduced  the  total  computing  time  at  most  20X  whereas  for  other  meshes  the  same 
strategy  actually  increased  the  computing  time.  It  was  decided  to  retain  the  more 
simple  algorithm  (see  Section  2)  in  whidn  K  is  updated  at  the  beoiming  of  each  load 
step  (the  computinig  times  shown  In  Table  2  are  for  that  algorithm). 

For  both  materials  60%  of  the  total  computing  time  was  involved  in  equation 
solving.  In  our  program,  the  banded  %jstem  of  linear  equations  is  solved  by  vectorized 
versions  of  the  UNPAdK  routines  and  SPBFA,  available  through  the  CRAY 

SSCILIB  library. 

5.  CONCLUSIONS. 

An  elastic-plastic  finite  element  method  which  uses  only  constant  strain 
triangular  elements  has  been  developed  for  a  problem  which  eimibits  a  stress 
singularity.  The  method  yields  a  complete  solution  for  stress,  strain,  and 
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displacement  which  is  accurate  even  at  points  which  are  very  close  to  the  singularity. 
The  method  does  not  require  any  9 priori  knowledge  of  the  form  of  the  sin^larity. 
The  load  is  applied  in  a  series  of  small  increments.  The  size  of  each  increment  Is 
determined  by  the  load  at  which  the  next  element  (out  of  all  those  that  are  currently 
elastic)  is  predicted  to  reach  yield.  The  nonlinear  material  behavior  is  accounted  for 
by  a  modified  Newton-Raphson  iterative  process  which  ensures  that  the  equilibrium 
equations  are  satisfied  at  the  end  of  each  load  increment.  The  solution  for  the 
particular  problem  of  a  centrally  cracked  plate  under  tensile  loading  in  a  state  of 
plane-stress  has  been  presented. 

The  finite  element  solution  was  obtained  by  using  a  mesh  in  which  the  size  of  the 
elements  decreases  in  a  geometric  series  as  the  crack  tip  is  approached  (see  Sec.  2). 
In  section  A,  it  was  shown  that  the  effect  of  different  element  arrangements  on  the 
elastic-plastic  solution  is  similar  to  that  observed  for  the  purely  elastic  solution 
(Ref.lH).  Features  are  included  in  the  finite  element  algorithm  which  significantly 
reduce  the  total  c.p.u.  time.  In  addition,  the  finite  elemeru  code  was  written  to  take 
advantage  of  the  vectorizing  capabilities  of  the  CRAY-1  computer,  thereby 
substantially  decreasing  the  computational  cost.  A  mesh  was  chosen,  for  our 
particular  problem,  wnich  provided  the  best  balance  between  accuracy  and 
computational  cost. 

Confidence  in  the  accuracy  of  the  solution  was  gained  from  a  comparison  with  the 
analytic  HRR  (Hutchinson-Rice-Rosengren)  [10,1 1,  and  121  solution  in  the  neighborhood 
of  the  crack  tip  under  conditions  of  small-scale  yielding.  It  was  shown  that  the  HRR 
solution  represents  the  behavior  of  the  crack  tip  stress-strain  field  provided  the 
applied  average  stress  does  not  exceed  about  one  half  of  the  yield  stress.  Under 
conditions  of  small-scale  yielding,  the  HRR  solution  characterized  the  stress-strain 
field  in  the  vicinity  of  the  crack  tip  only  over  a  very  small  portion  of  the  plastic  zone 
located  at  the  aaa  tip. 

The  method  presented  here  could  easily  be  extended  to  materials  which  obey  a 
piecewise  linear  stress-strain  law.  Yield  conditions  other  than  Von  Mises  could  also 
be  considered.  Other  fracture  problems  such  as  edge  cracks,  non-uniform  loading, 
shear  loading,  aacks  in  bending,  etc.  could  all  be  trivially  handled  by  changing  the 
boundary  conditions.  For  these  problems,  the  influence  of  geometry  and  boundary 
conditions  on  the  HRR  solution  could  be  studied.  This  has  applications  in  determining 
minimum  size  requirements  for  specimens  used  to  establish  a  "one  parameter”  fracture 
criterion  based  on  the  J-Inte^al.  The  method  could  be  modified  to  deal  with  cracks  in 
inhomogeneous  and/or  anisotropic  materials  such  as  composites. 

The  method  could  also  be  extended,  admittedly  not  without  some  effort,  to  study 
the  state  of  stress  at  the  tip  of  a  growing  crack  for  which  the  form  of  the  singularity 
is  not  well  understood  in  most  materials.  In  a  broader  context,  problems  involving 
other  forms  of  singularity  such  as  point  loads  and  reentrant  corners  could  also  be 
considered. 

In  this  paper,  the  emphasis  has  been  on  showing  that  a  finite  element  method  which 
uses  only  constant  strain  elements  can  provide  a  complete  elastic-plastic  solution  for 
a  state  of  plane-stress  even  in  the  region  of  high  stress  gradient  close  to  the 
singularity.  However  for  real  plates  the  state  of  stress  in  the  vicinity  of  the  crack 
tip  IS  fully  three  dimensional.  The  plane-stress  solution  can  be  expected  to  hold  only 
at  distances  from  the  aack  tip  which  exceed  the  plate  thickness.  This  means,  for  the 
particular  geometry  which  we  have  considered,  that  the  solution  can  have  physical 
meaning  only  over  a  range  of  about  r/a  >0.05  .  As  discussed  in  Section  3(b).  accurate 


solutions  can  be  obtained  over  this  range  using  a  coarse  mesh  (about  600  elements)  In 
relatively  small  c.  p.  u.  times.  It  has  been  sh^n  (Sec.  3(e))  that  large  strains  and 
rotations  are  predicted  in  a  neighbor^iood  surrounding  the  crack  tip.  the  size  of  which 
increases  as  the  allied  stress  is  increased.  The  assumption  of  small  strains  over  the 
range  r/a  >  0.05  holds  only  when  the  applied  stress  is  less  than  half  the  yield  stress. 

This  suggests  that  a  three  dimensional  finite  element  analysis  with  a  formulation 
to  account  Tor  large  strains  and  rotations  is  required  to  obtain  a  better  understanding 
of  the  fracture  process  in  the  region  close  to  the  crack  tip.  The  nonlinear  finite 
element  algorithm  and  mesh  ide%  Introduced  here  could  be  extended  to  the  three 
dimensional  problem  using  constant  strain  tetrahedral  elements.  The  increased  speed 
and  memory  capacity  of  Che  more  recent  models  of  the  CRAY  supercomputer  make  it 
reasonable  to  expect  that  a  full  three  dimensional  eiastlc-plastic  solution  can  be 
achieved  for  this  problem. 
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Table  II.  1.  Mesh  parameters. 
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Figure  1 1. 9.  Comparison  of  HRR  and  finite  element  stresses  along 
the  49®  ray. 
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